knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)
library(tidybayes)
library(ggplot2)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE 

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_final_with_plan_and_barriers_INTERACTION', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan', # todo: do we need this??
    'barriers_self_cw', 
    'facilitators_self_cw', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'barriers_self_cb', 
    'facilitators_self_cb', 
    'weartime_self_cb',
    
    # Interactions
    
    'persuasion_self_cw:barriers_self_cw',
    'barriers_self_cw:pressure_self_cw',
    'barriers_self_cw:pushing_self_cw',
    
    'persuasion_self_cw:facilitators_self_cw',
    'pressure_self_cw:facilitators_self_cw',
    'pushing_self_cw:facilitators_self_cw'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  
  'sd(barriers_self_cw)', 
  'sd(facilitators_self_cw)',
  
  'sd(persuasion_self_cw:barriers_self_cw)',
  'sd(barriers_self_cw:pressure_self_cw)',
  'sd(barriers_self_cw:pushing_self_cw)',
  
  'sd(persuasion_self_cw:facilitators_self_cw)',
  'sd(pressure_self_cw:facilitators_self_cw)',
  'sd(pushing_self_cw:facilitators_self_cw)',
  
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily individual's experienced persuasion",  
    "Daily partner's experienced persuasion", 
    "Daily individual's experienced pressure", 
    "Daily partner's experienced pressure", 
    "Daily individual's experienced pushing", 
    "Daily partner's experienced pushing", 
    "Day", 
    "Own actionplan",
    'Partner actionplan',
    "Daily barriers",
    "Daily facilitators",
    "Daily wear time",
    

    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean individual's experienced persuasion", 
    "Mean partner's experienced persuasion", 
    "Mean individual's experienced pressure", 
    "Mean partner's experienced pressure", 
    "Mean individual's experienced pushing", 
    "Mean partner's experienced pushing", 
    'Mean barriers',
    "Mean facilitators",
    "Mean wear time",
    
    # Interactions
    'persuasion_self_cw:barriers_self_cw',
    'barriers_self_cw:pressure_self_cw',
    'barriers_self_cw:pushing_self_cw',
    
    'persuasion_self_cw:facilitators_self_cw',
    'pressure_self_cw:facilitators_self_cw',
    'pushing_self_cw:facilitators_self_cw'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily individual's experienced persuasion)", 
  "sd(Daily partner's experienced persuasion)", # OR partner received
  "sd(Daily individual's experienced pressure)", 
  "sd(Daily partner's experienced pressure)", 
  "sd(Daily individual's experienced pushing)", 
  "sd(Daily partner's experienced pushing)", 
  
  'sd(barriers_self_cw)', 
  'sd(facilitators_self_cw)',
  
  'sd(persuasion_self_cw:barriers_self_cw)',
  'sd(barriers_self_cw:pressure_self_cw)',
  'sd(barriers_self_cw:pushing_self_cw)',
  
  'sd(persuasion_self_cw:facilitators_self_cw)',
  'sd(pressure_self_cw:facilitators_self_cw)',
  'sd(pushing_self_cw:facilitators_self_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,13),
  "Between-Person Effects" = c(14,22),
  "Interactions" = c(23, 28),
  "Random Effects" = c(29, 43), 
  "Additional Parameters" = c(44,44)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,13+5),
  "Between-Person Effects" = c(14+5,22+5),
  "Interactions" = c(23+5, 28+5),
  "Random Effects" = c(29+5, 43+5), 
  "Additional Parameters" = c(44+5,44+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day +  
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | dd | coupleID),
  
  hu = ~ persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | dd | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.30 [1.29, 1.30]         1.14      0.77
##       barriers_self_cw 1.15 [1.14, 1.15]         1.07      0.87
##       pressure_self_cw 1.08 [1.08, 1.09]         1.04      0.92
##        pushing_self_cw 1.09 [1.09, 1.09]         1.04      0.92
##   facilitators_self_cw 1.12 [1.12, 1.13]         1.06      0.89
##  persuasion_partner_cw 1.11 [1.10, 1.11]         1.05      0.90
##    pressure_partner_cw 1.05 [1.05, 1.05]         1.02      0.95
##     pushing_partner_cw 1.04 [1.04, 1.05]         1.02      0.96
##     persuasion_self_cb 1.05 [1.05, 1.05]         1.02      0.95
##  persuasion_partner_cb 3.68 [3.66, 3.70]         1.92      0.27
##       pressure_self_cb 3.60 [3.58, 3.62]         1.90      0.28
##    pressure_partner_cb 2.20 [2.19, 2.21]         1.48      0.45
##        pushing_self_cb 2.15 [2.14, 2.16]         1.47      0.47
##     pushing_partner_cb 3.01 [2.99, 3.02]         1.73      0.33
##       barriers_self_cb 3.05 [3.03, 3.07]         1.75      0.33
##   facilitators_self_cb 1.13 [1.12, 1.13]         1.06      0.89
##                    day 1.07 [1.07, 1.08]         1.04      0.93
##  Tolerance 95% CI
##      [0.77, 0.77]
##      [0.87, 0.88]
##      [0.92, 0.93]
##      [0.91, 0.92]
##      [0.89, 0.89]
##      [0.90, 0.91]
##      [0.95, 0.96]
##      [0.95, 0.96]
##      [0.95, 0.96]
##      [0.27, 0.27]
##      [0.28, 0.28]
##      [0.45, 0.46]
##      [0.46, 0.47]
##      [0.33, 0.33]
##      [0.33, 0.33]
##      [0.88, 0.89]
##      [0.93, 0.94]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         94%
##     lognormal          6%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         84%
##               beta-binomial          9%
##                   lognormal          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 47 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 6, observations = 1776, p-value =
## 0.14
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.003378378
## sample estimates:
## outlier frequency (expected: 0.00164414414414414 ) 
##                                        0.003378378
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report,
## rope_range = rope_range_continuous, : Coefficients were exponentiated.
## Double check if this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 3.41*** 1.09 [ 1.79, 6.47] 1.000 [0.84, 1.20] 0.001 55.719 Very Strong Evidence 1.000 25963 28327 50.01*** 4.20 [42.13, 59.09] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 8663 15562
Within-Person Effects
Daily individual’s experienced persuasion 1.87*** 0.28 [ 1.41, 2.61] 1.000 [0.84, 1.20] 0.001 >100 Overwhelming Evidence 1.000 32225 28628 1.00 0.03 [ 0.95, 1.06] 0.569 [0.93, 1.08] 0.993 0.010 Very Strong Evidence for Null 1.000 26933 31420
Daily partner’s experienced persuasion 1.40** 0.17 [ 1.12, 1.83] 0.998 [0.84, 1.20] 0.085 5.292 Moderate Evidence 1.000 32648 26113 1.02 0.02 [ 0.98, 1.07] 0.827 [0.93, 1.08] 0.992 0.013 Very Strong Evidence for Null 1.000 35756 32079
Daily individual’s experienced pressure 0.86 0.33 [ 0.39, 2.09] 0.653 [0.84, 1.20] 0.328 0.213 Moderate Evidence for Null 1.000 31627 24535 0.91 0.06 [ 0.78, 1.05] 0.911 [0.93, 1.08] 0.376 0.028 Very Strong Evidence for Null 1.000 30743 28087
Daily partner’s experienced pressure 1.69 0.72 [ 0.85, 5.72] 0.925 [0.84, 1.20] 0.162 0.451 Weak Evidence for Null 1.000 23866 22260 0.97 0.05 [ 0.88, 1.06] 0.757 [0.93, 1.08] 0.814 0.010 Very Strong Evidence for Null 1.000 40054 34187
Daily individual’s experienced pushing 1.24 0.28 [ 0.77, 2.01] 0.828 [0.84, 1.20] 0.389 0.188 Moderate Evidence for Null 1.000 29899 28023 0.99 0.03 [ 0.93, 1.06] 0.599 [0.93, 1.08] 0.972 0.009 Very Strong Evidence for Null 1.000 35123 35962
Daily partner’s experienced pushing 1.60** 0.30 [ 1.13, 2.54] 0.996 [0.84, 1.20] 0.048 3.367 Moderate Evidence 1.000 28148 22210 0.96 0.03 [ 0.90, 1.02] 0.911 [0.93, 1.08] 0.895 0.020 Very Strong Evidence for Null 1.000 37908 32301
Day 0.91 0.24 [ 0.55, 1.51] 0.636 [0.84, 1.20] 0.480 0.140 Moderate Evidence for Null 1.000 57086 30796 0.95 0.07 [ 0.83, 1.09] 0.745 [0.93, 1.08] 0.630 0.010 Very Strong Evidence for Null 1.000 54096 31338
Own actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.14 0.33 [ 0.64, 2.04] 0.670 [0.84, 1.20] 0.434 0.160 Moderate Evidence for Null 1.000 25852 28897 1.26*** 0.07 [ 1.13, 1.41] 1.000 [0.93, 1.08] 0.003 20.432 Strong Evidence 1.000 22243 28376
Daily facilitators 1.17 0.31 [ 0.68, 1.96] 0.724 [0.84, 1.20] 0.433 0.158 Moderate Evidence for Null 1.000 26886 28030 1.19** 0.07 [ 1.06, 1.33] 0.997 [0.93, 1.08] 0.054 0.970 Weak Evidence for Null 1.000 17913 23999
Daily wear time NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.07 0.73 [ 0.27, 4.02] 0.538 [0.84, 1.20] 0.207 0.338 Weak Evidence for Null 1.000 18162 26070 1.06 0.18 [ 0.75, 1.48] 0.633 [0.93, 1.08] 0.333 0.017 Very Strong Evidence for Null 1.000 16301 24978
Mean partner’s experienced persuasion 1.43 0.95 [ 0.38, 5.24] 0.706 [0.84, 1.20] 0.183 0.388 Weak Evidence for Null 1.000 18874 24215 0.96 0.16 [ 0.69, 1.32] 0.608 [0.93, 1.08] 0.347 0.016 Very Strong Evidence for Null 1.000 15922 24197
Mean individual’s experienced pressure 0.18 0.17 [ 0.03, 1.11] 0.969 [0.84, 1.20] 0.028 2.555 Weak Evidence 1.000 27633 30417 0.68 0.28 [ 0.29, 1.58] 0.823 [0.93, 1.08] 0.091 0.039 Strong Evidence for Null 1.001 6604 11985
Mean partner’s experienced pressure 0.19 0.18 [ 0.03, 1.17] 0.964 [0.84, 1.20] 0.031 2.363 Weak Evidence 1.000 25122 29791 0.49 0.21 [ 0.21, 1.15] 0.950 [0.93, 1.08] 0.035 0.102 Moderate Evidence for Null 1.001 6296 10725
Mean individual’s experienced pushing 1.32 1.26 [ 0.20, 8.70] 0.613 [0.84, 1.20] 0.143 0.485 Weak Evidence for Null 1.000 25432 29055 1.53 0.46 [ 0.84, 2.80] 0.920 [0.93, 1.08] 0.075 0.047 Strong Evidence for Null 1.000 8751 16967
Mean partner’s experienced pushing 1.45 1.41 [ 0.22, 9.51] 0.650 [0.84, 1.20] 0.137 0.517 Weak Evidence for Null 1.000 25870 29954 1.65 0.50 [ 0.91, 3.01] 0.950 [0.93, 1.08] 0.052 0.066 Strong Evidence for Null 1.000 8285 15361
Mean barriers 0.50** 0.13 [ 0.30, 0.83] 0.997 [0.84, 1.20] 0.023 5.207 Moderate Evidence 1.000 56783 30659 0.97 0.07 [ 0.83, 1.12] 0.674 [0.93, 1.08] 0.646 0.011 Very Strong Evidence for Null 1.000 40862 32230
Mean facilitators 1.24 0.23 [ 0.86, 1.80] 0.872 [0.84, 1.20] 0.415 0.175 Moderate Evidence for Null 1.000 54386 31332 1.08 0.05 [ 0.98, 1.19] 0.948 [0.93, 1.08] 0.479 0.049 Strong Evidence for Null 1.000 27594 29862
Mean wear time NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Interactions
persuasion_self_cw:barriers_self_cw 1.02 0.22 [ 0.65, 1.63] 0.528 [0.84, 1.20] 0.581 0.111 Moderate Evidence for Null 1.000 32162 28647 1.03 0.04 [ 0.96, 1.11] 0.777 [0.93, 1.08] 0.903 0.012 Very Strong Evidence for Null 1.000 38494 32777
barriers_self_cw:pressure_self_cw 0.34 0.22 [ 0.09, 1.21] 0.953 [0.84, 1.20] 0.054 1.436 Weak Evidence 1.000 34267 31019 0.95 0.11 [ 0.75, 1.20] 0.684 [0.93, 1.08] 0.441 0.012 Very Strong Evidence for Null 1.000 32615 31479
barriers_self_cw:pushing_self_cw 2.50* 1.01 [ 1.14, 6.03] 0.988 [0.84, 1.20] 0.029 2.808 Weak Evidence 1.000 26846 26558 0.94 0.05 [ 0.85, 1.05] 0.850 [0.93, 1.08] 0.632 0.018 Very Strong Evidence for Null 1.000 30501 30647
persuasion_self_cw:facilitators_self_cw 1.29 0.27 [ 0.87, 2.10] 0.890 [0.84, 1.20] 0.351 0.215 Moderate Evidence for Null 1.000 28751 26602 1.02 0.03 [ 0.96, 1.09] 0.735 [0.93, 1.08] 0.957 0.012 Very Strong Evidence for Null 1.000 27793 27637
pressure_self_cw:facilitators_self_cw 0.88 0.67 [ 0.18, 4.87] 0.568 [0.84, 1.20] 0.185 0.376 Weak Evidence for Null 1.000 29120 26540 0.94 0.10 [ 0.73, 1.18] 0.711 [0.93, 1.08] 0.462 0.014 Very Strong Evidence for Null 1.000 23694 20075
pushing_self_cw:facilitators_self_cw 1.54 0.57 [ 0.77, 3.92] 0.898 [0.84, 1.20] 0.190 0.383 Weak Evidence for Null 1.000 22894 22759 0.95 0.04 [ 0.88, 1.04] 0.879 [0.93, 1.08] 0.757 0.019 Very Strong Evidence for Null 1.000 32770 28974
Random Effects
sd(Intercept) 1.64 0.25 [1.23, 2.23] NA NA NA NA NA 1.000 23312 27758 0.32 0.05 [0.25, 0.43] NA NA NA NA NA 1.000 14324 22784
sd(Daily individual’s experienced persuasion) 0.21 0.17 [0.01, 0.61] NA NA NA NA NA 1.000 15715 20575 0.08 0.03 [0.03, 0.14] NA NA NA NA NA 1.000 9028 6570
sd(Daily partner’s experienced persuasion) 0.30 0.17 [0.03, 0.66] NA NA NA NA NA 1.000 12444 13969 0.06 0.02 [0.01, 0.11] NA NA NA NA NA 1.000 12117 8143
sd(Daily individual’s experienced pressure) 0.54 0.49 [0.02, 2.15] NA NA NA NA NA 1.000 17156 20291 0.08 0.08 [0.00, 0.33] NA NA NA NA NA 1.001 16273 21030
sd(Daily partner’s experienced pressure) 0.77 0.59 [0.04, 2.37] NA NA NA NA NA 1.000 13784 17928 0.05 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 20790 19556
sd(Daily individual’s experienced pushing) 0.70 0.36 [0.07, 1.53] NA NA NA NA NA 1.000 9718 10941 0.05 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 13373 17490
sd(Daily partner’s experienced pushing) 0.44 0.31 [0.03, 1.23] NA NA NA NA NA 1.000 11430 16990 0.06 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 11705 14178
sd(barriers_self_cw) 1.51 0.24 [1.10, 2.08] NA NA NA NA NA 1.000 30311 32165 0.23 0.05 [0.14, 0.35] NA NA NA NA NA 1.000 14821 17980
sd(facilitators_self_cw) 1.30 0.27 [0.84, 1.95] NA NA NA NA NA 1.000 24583 30170 0.27 0.05 [0.19, 0.39] NA NA NA NA NA 1.000 18417 26891
sd(persuasion_self_cw:barriers_self_cw) 0.58 0.33 [0.05, 1.26] NA NA NA NA NA 1.000 13106 13548 0.07 0.05 [0.00, 0.18] NA NA NA NA NA 1.001 9437 17501
sd(barriers_self_cw:pressure_self_cw) 0.67 0.61 [0.03, 2.58] NA NA NA NA NA 1.000 20813 20741 0.11 0.10 [0.01, 0.42] NA NA NA NA NA 1.000 21709 23743
sd(barriers_self_cw:pushing_self_cw) 1.26 0.49 [0.38, 2.39] NA NA NA NA NA 1.000 14916 14758 0.15 0.05 [0.05, 0.26] NA NA NA NA NA 1.000 14581 10756
sd(persuasion_self_cw:facilitators_self_cw) 0.59 0.27 [0.11, 1.23] NA NA NA NA NA 1.000 13497 10734 0.07 0.04 [0.00, 0.16] NA NA NA NA NA 1.000 9717 13565
sd(pressure_self_cw:facilitators_self_cw) 1.60 1.16 [0.10, 4.28] NA NA NA NA NA 1.000 13350 18455 0.12 0.11 [0.01, 0.51] NA NA NA NA NA 1.000 16895 22066
sd(pushing_self_cw:facilitators_self_cw) 0.94 0.62 [0.08, 2.45] NA NA NA NA NA 1.001 9901 14021 0.05 0.04 [0.00, 0.17] NA NA NA NA NA 1.000 19197 21830
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.62 0.01 [0.59, 0.65] NA NA NA NA NA 1.000 37840 29108
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.77), b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.75), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.83), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.71). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Hurdle part of the model on the left, non-zero part towards the right side of the table

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Scanning ttf files in C:\Windows/Fonts, C:\Users\pascku\AppData\Local/Microsoft/Windows/Fonts ...
## Extracting .afm files from .ttf files...
## C:\Windows\Fonts\Candara.ttf : Candara already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarab.ttf : Candara-Bold already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarai.ttf : Candara-Italic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaral.ttf : Candara-Light already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarali.ttf : Candara-LightItalic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaraz.ttf : Candara-BoldItalic already registered in fonts database. Skipping.
## Found FontName for 0 fonts.
## Scanning afm files in C:/Users/pascku/AppData/Local/R/cache/R/renv/cache/v5/windows/R-4.4/x86_64-w64-mingw32/extrafontdb/1.0/a861555ddec7451c653b40e713166c6f/extrafontdb/metrics
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 98 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 191 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.013

$persuasion_partner_cw

## Warning: Removed 118 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 240 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.00825

$pressure_self_cw

## Warning: Removed 76 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 117 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.0185

$pressure_partner_cw

## Warning: Removed 75 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 100 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.0326

$pushing_self_cw

## Warning: Removed 100 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 58 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.0136

$pushing_partner_cw

## Warning: Removed 35 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 75 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.014

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Comparing effect size of pressure and pushing

hypothesis(pa_sub, "pressure_self_cw < pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... < 0    -0.09      0.09    -0.23     0.05       5.71
##   Post.Prob Star
## 1      0.85     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  control = list(adapt_delta = 0.99),
  iter = iterations + 200,
  warmup = warmup + 200,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.20 [1.16, 1.24]         1.09      0.84
##       barriers_self_cw 1.20 [1.17, 1.24]         1.10      0.83
##       pressure_self_cw 1.17 [1.13, 1.21]         1.08      0.86
##        pushing_self_cw 1.26 [1.22, 1.30]         1.12      0.80
##   facilitators_self_cw 1.17 [1.14, 1.21]         1.08      0.85
##  persuasion_partner_cw 1.14 [1.11, 1.18]         1.07      0.88
##    pressure_partner_cw 1.10 [1.08, 1.14]         1.05      0.91
##     pushing_partner_cw 1.16 [1.12, 1.20]         1.08      0.86
##       pressure_self_cb 3.80 [3.62, 3.98]         1.95      0.26
##    pressure_partner_cb 4.01 [3.83, 4.21]         2.00      0.25
##       barriers_self_cb 1.18 [1.15, 1.22]         1.09      0.85
##   facilitators_self_cb 1.16 [1.13, 1.20]         1.08      0.86
##                    day 1.06 [1.04, 1.10]         1.03      0.94
##       weartime_self_cw 1.06 [1.04, 1.10]         1.03      0.94
##       weartime_self_cb 1.21 [1.18, 1.26]         1.10      0.82
##  Tolerance 95% CI
##      [0.81, 0.86]
##      [0.80, 0.86]
##      [0.83, 0.88]
##      [0.77, 0.82]
##      [0.83, 0.88]
##      [0.85, 0.90]
##      [0.88, 0.93]
##      [0.84, 0.89]
##      [0.25, 0.28]
##      [0.24, 0.26]
##      [0.82, 0.87]
##      [0.83, 0.89]
##      [0.91, 0.96]
##      [0.91, 0.96]
##      [0.80, 0.85]
## 
## Moderate Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 7.76 [7.38, 8.17]         2.79      0.13
##  persuasion_partner_cb 8.17 [7.77, 8.60]         2.86      0.12
##        pushing_self_cb 5.88 [5.60, 6.18]         2.42      0.17
##     pushing_partner_cb 5.74 [5.47, 6.03]         2.40      0.17
##  Tolerance 95% CI
##      [0.12, 0.14]
##      [0.12, 0.13]
##      [0.16, 0.18]
##      [0.17, 0.18]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 11 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 14, observations = 1594, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0002979925 0.0056461731
## sample estimates:
## outlier frequency (expected: 0.003099121706399 ) 
##                                      0.008782936
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 123.06*** 7.77 [108.16, 139.27] 1.000 [0.94, 1.06] 0.000 >100 Overwhelming Evidence 1.000 8031 14512
Within-Person Effects
Daily individual’s experienced persuasion 1.02 0.02 [ 0.98, 1.06] 0.864 [0.94, 1.06] 0.985 0.013 Very Strong Evidence for Null 1.000 34459 31842
Daily partner’s experienced persuasion 1.03 0.02 [ 0.99, 1.07] 0.931 [0.94, 1.06] 0.969 0.021 Very Strong Evidence for Null 1.000 38113 30960
Daily individual’s experienced pressure 0.97 0.05 [ 0.87, 1.07] 0.759 [0.94, 1.06] 0.672 0.009 Very Strong Evidence for Null 1.000 40938 32241
Daily partner’s experienced pressure 0.97 0.04 [ 0.90, 1.05] 0.793 [0.94, 1.06] 0.797 0.009 Very Strong Evidence for Null 1.000 48631 31224
Daily individual’s experienced pushing 1.02 0.03 [ 0.96, 1.08] 0.734 [0.94, 1.06] 0.939 0.009 Very Strong Evidence for Null 1.000 37349 33850
Daily partner’s experienced pushing 1.00 0.02 [ 0.96, 1.05] 0.557 [0.94, 1.06] 0.988 0.006 Very Strong Evidence for Null 1.000 43220 32281
Day 0.96 0.05 [ 0.87, 1.06] 0.771 [0.94, 1.06] 0.666 0.008 Very Strong Evidence for Null 1.000 60635 30163
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.01 0.03 [ 0.94, 1.07] 0.562 [0.94, 1.06] 0.939 0.008 Very Strong Evidence for Null 1.000 28586 29795
Daily facilitators 1.04 0.03 [ 0.99, 1.10] 0.945 [0.94, 1.06] 0.778 0.031 Strong Evidence for Null 1.000 31594 27535
Daily wear time 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.06] 1.000 2.621 Weak Evidence 1.000 66527 31035
Between-Person Effects
Mean individual’s experienced persuasion 0.99 0.17 [ 0.69, 1.40] 0.533 [0.94, 1.06] 0.274 0.016 Very Strong Evidence for Null 1.001 6637 13062
Mean partner’s experienced persuasion 0.86 0.15 [ 0.61, 1.23] 0.796 [0.94, 1.06] 0.196 0.023 Very Strong Evidence for Null 1.001 6727 12797
Mean individual’s experienced pressure 1.12 0.21 [ 0.77, 1.64] 0.730 [0.94, 1.06] 0.216 0.012 Very Strong Evidence for Null 1.001 9668 17120
Mean partner’s experienced pressure 1.04 0.19 [ 0.72, 1.49] 0.579 [0.94, 1.06] 0.261 0.011 Very Strong Evidence for Null 1.001 8766 16544
Mean individual’s experienced pushing 1.11 0.27 [ 0.69, 1.82] 0.668 [0.94, 1.06] 0.186 0.014 Very Strong Evidence for Null 1.001 9898 17862
Mean partner’s experienced pushing 1.27 0.31 [ 0.78, 2.05] 0.836 [0.94, 1.06] 0.127 0.022 Very Strong Evidence for Null 1.001 9438 16390
Mean barriers 0.80*** 0.04 [ 0.72, 0.89] 1.000 [0.94, 1.06] 0.001 26.021 Strong Evidence 1.000 52877 32143
Mean facilitators 0.96 0.04 [ 0.89, 1.03] 0.873 [0.94, 1.06] 0.704 0.019 Very Strong Evidence for Null 1.000 34078 32619
Mean wear time 1.00* 0.00 [ 1.00, 1.00] 0.988 [0.94, 1.06] 1.000 0.124 Moderate Evidence for Null 1.000 31806 31751
Interactions
persuasion_self_cw:barriers_self_cw 0.97 0.02 [ 0.93, 1.02] 0.882 [0.94, 1.06] 0.905 0.013 Very Strong Evidence for Null 1.000 53504 32603
barriers_self_cw:pressure_self_cw 0.96 0.08 [ 0.81, 1.14] 0.684 [0.94, 1.06] 0.495 0.008 Very Strong Evidence for Null 1.000 34654 33293
barriers_self_cw:pushing_self_cw 1.02 0.03 [ 0.95, 1.09] 0.696 [0.94, 1.06] 0.894 0.007 Very Strong Evidence for Null 1.000 43986 31833
persuasion_self_cw:facilitators_self_cw 0.99 0.02 [ 0.95, 1.04] 0.629 [0.94, 1.06] 0.994 0.008 Very Strong Evidence for Null 1.000 38120 29806
pressure_self_cw:facilitators_self_cw 0.92 0.10 [ 0.69, 1.14] 0.798 [0.94, 1.06] 0.328 0.016 Very Strong Evidence for Null 1.000 19016 20911
pushing_self_cw:facilitators_self_cw 1.03 0.03 [ 0.96, 1.09] 0.793 [0.94, 1.06] 0.868 0.011 Very Strong Evidence for Null 1.000 36408 28457
Random Effects
sd(Intercept) 0.34 0.05 [0.26, 0.45] NA NA NA NA NA 1.000 10470 18026
sd(Daily individual’s experienced persuasion) 0.05 0.02 [0.01, 0.09] NA NA NA NA NA 1.000 12921 10061
sd(Daily partner’s experienced persuasion) 0.05 0.02 [0.01, 0.10] NA NA NA NA NA 1.000 10508 9922
sd(Daily individual’s experienced pressure) 0.05 0.05 [0.00, 0.20] NA NA NA NA NA 1.000 16988 19408
sd(Daily partner’s experienced pressure) 0.03 0.03 [0.00, 0.12] NA NA NA NA NA 1.000 22902 19193
sd(Daily individual’s experienced pushing) 0.08 0.04 [0.01, 0.15] NA NA NA NA NA 1.000 10344 11953
sd(Daily partner’s experienced pushing) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 12602 15278
sd(barriers_self_cw) 0.13 0.03 [0.06, 0.20] NA NA NA NA NA 1.000 11395 8797
sd(facilitators_self_cw) 0.07 0.04 [0.01, 0.15] NA NA NA NA NA 1.000 9088 10195
sd(persuasion_self_cw:barriers_self_cw) 0.03 0.02 [0.00, 0.09] NA NA NA NA NA 1.000 19050 20948
sd(barriers_self_cw:pressure_self_cw) 0.07 0.07 [0.00, 0.31] NA NA NA NA NA 1.000 20796 21799
sd(barriers_self_cw:pushing_self_cw) 0.03 0.03 [0.00, 0.10] NA NA NA NA NA 1.000 21260 21607
sd(persuasion_self_cw:facilitators_self_cw) 0.03 0.02 [0.00, 0.09] NA NA NA NA NA 1.000 13572 18030
sd(pressure_self_cw:facilitators_self_cw) 0.24 0.16 [0.02, 0.62] NA NA NA NA NA 1.000 10644 14130
sd(pushing_self_cw:facilitators_self_cw) 0.04 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 15627 18450
Additional Parameters
sigma 0.53 0.01 [0.51, 0.55] NA NA NA NA NA 1.000 40490 31103
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.87), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.72), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.78). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.24 [1.20,  1.29]         1.11      0.81
##       barriers_self_cw 1.12 [1.09,  1.16]         1.06      0.89
##       pressure_self_cw 1.14 [1.11,  1.18]         1.07      0.87
##        pushing_self_cw 1.30 [1.26,  1.34]         1.14      0.77
##   facilitators_self_cw 1.09 [1.06,  1.13]         1.04      0.92
##  persuasion_partner_cw 1.14 [1.11,  1.19]         1.07      0.87
##    pressure_partner_cw 1.09 [1.06,  1.13]         1.04      0.92
##     pushing_partner_cw 1.13 [1.10,  1.17]         1.06      0.88
##       pressure_self_cb 5.00 [4.75,  5.25]         2.24      0.20
##    pressure_partner_cb 4.87 [4.63,  5.12]         2.21      0.21
##       barriers_self_cb 1.12 [1.09,  1.16]         1.06      0.89
##   facilitators_self_cb 1.09 [1.06,  1.13]         1.04      0.92
##                    day 1.08 [1.06,  1.12]         1.04      0.92
##  Tolerance 95% CI
##      [0.78, 0.83]
##      [0.86, 0.92]
##      [0.84, 0.90]
##      [0.74, 0.80]
##      [0.89, 0.94]
##      [0.84, 0.90]
##      [0.89, 0.94]
##      [0.85, 0.91]
##      [0.19, 0.21]
##      [0.20, 0.22]
##      [0.86, 0.92]
##      [0.88, 0.94]
##      [0.89, 0.95]
## 
## Moderate Correlation
## 
##                Term  VIF    VIF 95% CI Increased SE Tolerance
##     pushing_self_cb 7.80 [7.41,  8.22]         2.79      0.13
##  pushing_partner_cb 7.80 [7.41,  8.22]         2.79      0.13
##  Tolerance 95% CI
##      [0.12, 0.13]
##      [0.12, 0.13]
## 
## High Correlation
## 
##                   Term   VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 10.41 [9.88, 10.98]         3.23      0.10
##  persuasion_partner_cb 10.26 [9.73, 10.81]         3.20      0.10
##  Tolerance 95% CI
##      [0.09, 0.10]
##      [0.09, 0.10]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         41%
##        normal         41%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 10 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 9, observations = 1776, p-value =
## 0.01064
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.002319754 0.009597934
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.005067568
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.
## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF    VIF 95% CI Increased SE Tolerance
persuasion_self_cw 1.24 [1.20,  1.29]         1.11      0.81
  barriers_self_cw 1.12 [1.09,  1.16]         1.06      0.89
  pressure_self_cw 1.14 [1.11,  1.18]         1.07      0.87
   pushing_self_cw 1.30 [1.26,  1.34]         1.14      0.77

facilitators_self_cw 1.09 [1.06, 1.13] 1.04 0.92 persuasion_partner_cw 1.14 [1.11, 1.19] 1.07 0.87 pressure_partner_cw 1.09 [1.06, 1.13] 1.04 0.92 pushing_partner_cw 1.13 [1.10, 1.17] 1.06 0.88 pressure_self_cb 5.00 [4.75, 5.25] 2.24 0.20 pressure_partner_cb 4.87 [4.63, 5.12] 2.21 0.21 barriers_self_cb 1.12 [1.09, 1.16] 1.06 0.89 facilitators_self_cb 1.09 [1.06, 1.13] 1.04 0.92 day 1.08 [1.06, 1.12] 1.04 0.92 Tolerance 95% CI [0.78, 0.83] [0.86, 0.92] [0.84, 0.90] [0.74, 0.80] [0.89, 0.94] [0.84, 0.90] [0.89, 0.94] [0.85, 0.91] [0.19, 0.21] [0.20, 0.22] [0.86, 0.92] [0.88, 0.94] [0.89, 0.95]

Moderate Correlation

           Term  VIF    VIF 95% CI Increased SE Tolerance
pushing_self_cb 7.80 [7.41,  8.22]         2.79      0.13

pushing_partner_cb 7.80 [7.41, 8.22] 2.79 0.13 Tolerance 95% CI [0.12, 0.13] [0.12, 0.13]

High Correlation

              Term   VIF    VIF 95% CI Increased SE Tolerance
persuasion_self_cb 10.41 [9.88, 10.98]         3.23      0.10

persuasion_partner_cb 10.26 [9.73, 10.81] 3.20 0.10 Tolerance 95% CI [0.09, 0.10] [0.09, 0.10]

## Sampling priors, please wait...
summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.79*** 0.11 [ 3.58, 4.01] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 8034 15527
Within-Person Effects
Daily individual’s experienced persuasion 0.01 0.03 [-0.04, 0.06] 0.615 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 58341 33448
Daily partner’s experienced persuasion 0.01 0.02 [-0.04, 0.06] 0.663 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 58314 31804
Daily individual’s experienced pressure -0.02 0.07 [-0.16, 0.12] 0.619 [-0.11, 0.11] 0.859 0.005 Very Strong Evidence for Null 1.000 47243 34060
Daily partner’s experienced pressure -0.05 0.06 [-0.18, 0.07] 0.810 [-0.11, 0.11] 0.825 0.008 Very Strong Evidence for Null 1.000 44211 31698
Daily individual’s experienced pushing 0.02 0.04 [-0.06, 0.10] 0.705 [-0.11, 0.11] 0.980 0.006 Very Strong Evidence for Null 1.000 43183 35615
Daily partner’s experienced pushing 0.09* 0.04 [ 0.02, 0.16] 0.990 [-0.11, 0.11] 0.690 0.097 Strong Evidence for Null 1.000 47310 30460
Day 0.26*** 0.08 [ 0.11, 0.41] 1.000 [-0.11, 0.11] 0.025 1.134 Weak Evidence 1.000 67179 31237
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers -0.12 0.06 [-0.25, 0.01] 0.969 [-0.11, 0.11] 0.431 0.049 Strong Evidence for Null 1.000 20011 25278
Daily facilitators 0.11** 0.03 [ 0.05, 0.17] 0.999 [-0.11, 0.11] 0.502 1.022 Weak Evidence 1.000 55329 29865
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.13 0.30 [-0.47, 0.72] 0.667 [-0.11, 0.11] 0.265 0.015 Very Strong Evidence for Null 1.001 7278 14121
Mean partner’s experienced persuasion 0.42 0.29 [-0.18, 1.00] 0.918 [-0.11, 0.11] 0.112 0.040 Strong Evidence for Null 1.001 7123 13430
Mean individual’s experienced pressure -0.07 0.30 [-0.67, 0.53] 0.594 [-0.11, 0.11] 0.275 0.009 Very Strong Evidence for Null 1.000 8986 17404
Mean partner’s experienced pressure -0.32 0.30 [-0.91, 0.28] 0.854 [-0.11, 0.11] 0.166 0.016 Very Strong Evidence for Null 1.000 8783 16707
Mean individual’s experienced pushing 0.23 0.42 [-0.61, 1.08] 0.705 [-0.11, 0.11] 0.178 0.014 Very Strong Evidence for Null 1.001 10000 16141
Mean partner’s experienced pushing -0.06 0.42 [-0.90, 0.79] 0.557 [-0.11, 0.11] 0.208 0.012 Very Strong Evidence for Null 1.000 9929 16149
Mean barriers -0.50*** 0.08 [-0.65, -0.34] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 72710 31917
Mean facilitators 0.29*** 0.05 [ 0.18, 0.39] 1.000 [-0.11, 0.11] 0.001 >100 Overwhelming Evidence 1.000 43379 32312
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Interactions
persuasion_self_cw:barriers_self_cw 0.00 0.04 [-0.08, 0.07] 0.529 [-0.11, 0.11] 0.993 0.005 Very Strong Evidence for Null 1.000 50722 33188
barriers_self_cw:pressure_self_cw -0.14 0.16 [-0.47, 0.20] 0.808 [-0.11, 0.11] 0.364 0.011 Very Strong Evidence for Null 1.000 31974 31759
barriers_self_cw:pushing_self_cw 0.08 0.05 [-0.02, 0.18] 0.950 [-0.11, 0.11] 0.701 0.019 Very Strong Evidence for Null 1.000 55985 33379
persuasion_self_cw:facilitators_self_cw 0.00 0.03 [-0.06, 0.06] 0.530 [-0.11, 0.11] 0.999 0.005 Very Strong Evidence for Null 1.000 51081 31523
pressure_self_cw:facilitators_self_cw 0.00 0.13 [-0.24, 0.31] 0.502 [-0.11, 0.11] 0.611 0.007 Very Strong Evidence for Null 1.000 27612 25163
pushing_self_cw:facilitators_self_cw 0.02 0.05 [-0.08, 0.13] 0.673 [-0.11, 0.11] 0.935 0.007 Very Strong Evidence for Null 1.000 34592 31350
Random Effects
sd(Intercept) 0.62 0.08 [0.49, 0.81] NA NA NA NA NA 1.000 9965 19755
sd(Daily individual’s experienced persuasion) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 13775 19005
sd(Daily partner’s experienced persuasion) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 16581 19683
sd(Daily individual’s experienced pressure) 0.07 0.06 [0.00, 0.25] NA NA NA NA NA 1.000 23395 22857
sd(Daily partner’s experienced pressure) 0.09 0.07 [0.00, 0.29] NA NA NA NA NA 1.000 16353 19990
sd(Daily individual’s experienced pushing) 0.09 0.05 [0.01, 0.20] NA NA NA NA NA 1.001 11251 12691
sd(Daily partner’s experienced pushing) 0.09 0.05 [0.01, 0.19] NA NA NA NA NA 1.000 13536 16131
sd(barriers_self_cw) 0.32 0.05 [0.23, 0.43] NA NA NA NA NA 1.000 19603 27622
sd(facilitators_self_cw) 0.05 0.04 [0.00, 0.16] NA NA NA NA NA 1.000 13807 20577
sd(persuasion_self_cw:barriers_self_cw) 0.07 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 14309 19001
sd(barriers_self_cw:pressure_self_cw) 0.41 0.18 [0.08, 0.85] NA NA NA NA NA 1.000 14368 12763
sd(barriers_self_cw:pushing_self_cw) 0.05 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 20198 24464
sd(persuasion_self_cw:facilitators_self_cw) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 17608 20566
sd(pressure_self_cw:facilitators_self_cw) 0.16 0.13 [0.01, 0.55] NA NA NA NA NA 1.000 19179 20271
sd(pushing_self_cw:facilitators_self_cw) 0.11 0.06 [0.01, 0.26] NA NA NA NA NA 1.000 14431 15173
Additional Parameters
sigma 0.84 0.02 [0.82, 0.88] NA NA NA NA NA 1.000 53758 30528
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.73), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.73), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.75), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.85). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##       pressure_self_cw 4.08 [3.79,  4.40]         2.02      0.24
##        pushing_self_cw 1.37 [1.31,  1.45]         1.17      0.73
##   facilitators_self_cw 1.27 [1.21,  1.34]         1.13      0.79
##  persuasion_partner_cw 1.69 [1.60,  1.80]         1.30      0.59
##    pressure_partner_cw 1.12 [1.08,  1.19]         1.06      0.89
##     pushing_partner_cw 1.29 [1.23,  1.37]         1.14      0.77
##     persuasion_self_cb 1.91 [1.80,  2.04]         1.38      0.52
##  persuasion_partner_cb 1.22 [1.17,  1.29]         1.11      0.82
##       pressure_self_cb 1.06 [1.03,  1.13]         1.03      0.94
##    pressure_partner_cb 1.23 [1.18,  1.30]         1.11      0.81
##       barriers_self_cb 4.96 [4.60,  5.36]         2.23      0.20
##   facilitators_self_cb 4.96 [4.60,  5.36]         2.23      0.20
##                    day 2.14 [2.01,  2.28]         1.46      0.47
##  Tolerance 95% CI
##      [0.23, 0.26]
##      [0.69, 0.77]
##      [0.75, 0.83]
##      [0.55, 0.62]
##      [0.84, 0.93]
##      [0.73, 0.81]
##      [0.49, 0.56]
##      [0.77, 0.86]
##      [0.88, 0.97]
##      [0.77, 0.85]
##      [0.19, 0.22]
##      [0.19, 0.22]
##      [0.44, 0.50]
## 
## Moderate Correlation
## 
##                Term  VIF    VIF 95% CI Increased SE Tolerance
##  persuasion_self_cw 8.57 [7.91,  9.28]         2.93      0.12
##     pushing_self_cb 5.03 [4.66,  5.43]         2.24      0.20
##  pushing_partner_cb 6.21 [5.75,  6.72]         2.49      0.16
##  Tolerance 95% CI
##      [0.11, 0.13]
##      [0.18, 0.21]
##      [0.15, 0.17]
## 
## High Correlation
## 
##              Term   VIF    VIF 95% CI Increased SE Tolerance
##  barriers_self_cw 10.36 [9.56, 11.24]         3.22      0.10
##  Tolerance 95% CI
##      [0.09, 0.10]
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'print': Predictive errors are not defined for ordinal or categorical models.
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 36 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 486, p-value = 0.38
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002057613
## sample estimates:
## outlier frequency (expected: 0.000432098765432099 ) 
##                                         0.002057613
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.
## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF    VIF 95% CI Increased SE Tolerance
  pressure_self_cw 4.08 [3.79,  4.40]         2.02      0.24
   pushing_self_cw 1.37 [1.31,  1.45]         1.17      0.73

facilitators_self_cw 1.27 [1.21, 1.34] 1.13 0.79 persuasion_partner_cw 1.69 [1.60, 1.80] 1.30 0.59 pressure_partner_cw 1.12 [1.08, 1.19] 1.06 0.89 pushing_partner_cw 1.29 [1.23, 1.37] 1.14 0.77 persuasion_self_cb 1.91 [1.80, 2.04] 1.38 0.52 persuasion_partner_cb 1.22 [1.17, 1.29] 1.11 0.82 pressure_self_cb 1.06 [1.03, 1.13] 1.03 0.94 pressure_partner_cb 1.23 [1.18, 1.30] 1.11 0.81 barriers_self_cb 4.96 [4.60, 5.36] 2.23 0.20 facilitators_self_cb 4.96 [4.60, 5.36] 2.23 0.20 day 2.14 [2.01, 2.28] 1.46 0.47 Tolerance 95% CI [0.23, 0.26] [0.69, 0.77] [0.75, 0.83] [0.55, 0.62] [0.84, 0.93] [0.73, 0.81] [0.49, 0.56] [0.77, 0.86] [0.88, 0.97] [0.77, 0.85] [0.19, 0.22] [0.19, 0.22] [0.44, 0.50]

Moderate Correlation

           Term  VIF    VIF 95% CI Increased SE Tolerance

persuasion_self_cw 8.57 [7.91, 9.28] 2.93 0.12 pushing_self_cb 5.03 [4.66, 5.43] 2.24 0.20 pushing_partner_cb 6.21 [5.75, 6.72] 2.49 0.16 Tolerance 95% CI [0.11, 0.13] [0.18, 0.21] [0.15, 0.17]

High Correlation

         Term   VIF    VIF 95% CI Increased SE Tolerance

barriers_self_cw 10.36 [9.56, 11.24] 3.22 0.10 Tolerance 95% CI [0.09, 0.10]

## Sampling priors, please wait...
summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 5.24*** 2.33 [ 2.23, 13.19] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 33731 29944
Intercept[2] 13.11*** 6.13 [ 5.42, 34.96] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 31987 29285
Intercept[3] 43.61*** 22.10 [ 17.01, 128.11] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 29869 29288
Intercept[4] 254.83*** 155.72 [ 83.32, 917.86] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 27685 28142
Intercept[5] 13351.51*** 16831.35 [1628.97, 285902.89] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 39002 29466
Within-Person Effects
Daily individual’s experienced persuasion 0.64** 0.10 [ 0.46, 0.86] 0.999 [0.84, 1.20] 0.039 7.369 Moderate Evidence 1.000 32737 28343
Daily partner’s experienced persuasion 1.10 0.17 [ 0.80, 1.48] 0.741 [0.84, 1.20] 0.659 0.088 Strong Evidence for Null 1.000 36870 29745
Daily individual’s experienced pressure 2.21* 0.71 [ 1.08, 4.26] 0.983 [0.84, 1.20] 0.037 1.223 Weak Evidence 1.000 26836 26589
Daily partner’s experienced pressure 1.13 0.43 [ 0.41, 2.46] 0.628 [0.84, 1.20] 0.337 0.091 Strong Evidence for Null 1.000 26133 20054
Daily individual’s experienced pushing 1.30 0.20 [ 0.93, 1.78] 0.944 [0.84, 1.20] 0.302 0.271 Moderate Evidence for Null 1.000 38983 28672
Daily partner’s experienced pushing 0.92 0.16 [ 0.65, 1.32] 0.684 [0.84, 1.20] 0.642 0.075 Strong Evidence for Null 1.000 41776 30483
Day 1.33 0.76 [ 0.43, 4.11] 0.688 [0.84, 1.20] 0.217 0.071 Strong Evidence for Null 1.000 45608 32116
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.28 0.51 [ 0.58, 2.94] 0.737 [0.84, 1.20] 0.290 0.086 Strong Evidence for Null 1.000 43277 33473
Daily facilitators 1.15 0.51 [ 0.44, 2.67] 0.618 [0.84, 1.20] 0.294 0.093 Strong Evidence for Null 1.000 29957 29903
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.49 1.31 [ 0.25, 9.25] 0.674 [0.84, 1.20] 0.149 0.100 Strong Evidence for Null 1.000 22625 27377
Mean partner’s experienced persuasion 1.94 1.88 [ 0.31, 15.31] 0.762 [0.84, 1.20] 0.120 0.114 Moderate Evidence for Null 1.000 22074 25540
Mean individual’s experienced pressure 2.88 3.01 [ 0.37, 24.66] 0.848 [0.84, 1.20] 0.083 0.128 Moderate Evidence for Null 1.000 26998 30390
Mean partner’s experienced pressure 1.08 1.09 [ 0.13, 6.97] 0.531 [0.84, 1.20] 0.142 0.078 Strong Evidence for Null 1.000 21451 25485
Mean individual’s experienced pushing 1.02 1.29 [ 0.08, 14.46] 0.507 [0.84, 1.20] 0.115 0.086 Strong Evidence for Null 1.000 23211 27282
Mean partner’s experienced pushing 0.02** 0.03 [ 0.00, 0.33] 0.996 [0.84, 1.20] 0.003 3.686 Moderate Evidence 1.000 24552 26788
Mean barriers 0.89 0.63 [ 0.21, 3.48] 0.564 [0.84, 1.20] 0.200 0.081 Strong Evidence for Null 1.000 32718 28853
Mean facilitators 0.60 0.28 [ 0.23, 1.49] 0.871 [0.84, 1.20] 0.168 0.194 Moderate Evidence for Null 1.000 28304 29905
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Interactions
persuasion_self_cw:barriers_self_cw 0.69 0.15 [ 0.43, 1.04] 0.961 [0.84, 1.20] 0.182 0.365 Weak Evidence for Null 1.000 32999 29006
barriers_self_cw:pressure_self_cw 1.71 0.75 [ 0.70, 4.17] 0.888 [0.84, 1.20] 0.153 0.157 Moderate Evidence for Null 1.000 35588 33208
barriers_self_cw:pushing_self_cw 0.69 0.16 [ 0.42, 1.08] 0.947 [0.84, 1.20] 0.195 0.264 Moderate Evidence for Null 1.000 35035 29498
persuasion_self_cw:facilitators_self_cw 0.91 0.18 [ 0.62, 1.35] 0.682 [0.84, 1.20] 0.586 0.097 Strong Evidence for Null 1.000 33388 30290
pressure_self_cw:facilitators_self_cw 1.68 1.00 [ 0.55, 6.63] 0.820 [0.84, 1.20] 0.168 0.173 Moderate Evidence for Null 1.000 18589 20585
pushing_self_cw:facilitators_self_cw 0.65 0.16 [ 0.40, 1.12] 0.945 [0.84, 1.20] 0.150 0.431 Weak Evidence for Null 1.000 28963 23880
Random Effects
sd(Intercept) 1.30 0.34 [0.74, 2.12] NA NA NA NA NA 1.000 13995 22011
sd(Daily individual’s experienced persuasion) 0.34 0.23 [0.02, 0.81] NA NA NA NA NA 1.000 8114 15291
sd(Daily partner’s experienced persuasion) 0.24 0.18 [0.01, 0.64] NA NA NA NA NA 1.000 15429 19416
sd(Daily individual’s experienced pressure) 0.74 0.42 [0.07, 1.75] NA NA NA NA NA 1.000 12028 12872
sd(Daily partner’s experienced pressure) 0.75 0.61 [0.04, 2.34] NA NA NA NA NA 1.000 11524 18542
sd(Daily individual’s experienced pushing) 0.26 0.22 [0.01, 0.79] NA NA NA NA NA 1.000 11607 19536
sd(Daily partner’s experienced pushing) 0.20 0.18 [0.01, 0.75] NA NA NA NA NA 1.000 18835 21947
sd(barriers_self_cw) 0.75 0.52 [0.04, 1.86] NA NA NA NA NA 1.000 10204 16342
sd(facilitators_self_cw) 1.15 0.44 [0.23, 2.08] NA NA NA NA NA 1.001 9945 7825
sd(persuasion_self_cw:barriers_self_cw) 0.37 0.25 [0.02, 0.93] NA NA NA NA NA 1.000 13473 16706
sd(barriers_self_cw:pressure_self_cw) 0.41 0.37 [0.02, 1.62] NA NA NA NA NA 1.000 21785 20837
sd(barriers_self_cw:pushing_self_cw) 0.30 0.25 [0.01, 0.94] NA NA NA NA NA 1.000 14225 20193
sd(persuasion_self_cw:facilitators_self_cw) 0.21 0.18 [0.01, 0.70] NA NA NA NA NA 1.001 14719 19964
sd(pressure_self_cw:facilitators_self_cw) 1.19 0.62 [0.16, 2.69] NA NA NA NA NA 1.000 13974 13085
sd(pushing_self_cw:facilitators_self_cw) 0.35 0.29 [0.02, 1.12] NA NA NA NA NA 1.000 13347 21462
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2]
##   (r = 0.79), b_Intercept[4] and b_Intercept[3] (r = 0.85),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.76). This
##   might lead to inappropriate results. See 'Details' in '?rope'.

#conditional_spaghetti(
#  reactance_ordinal, 
#  effects = c('persuasion_self_cw', 'pressure_self_cw')
#  , group_var = 'coupleID'
#  #, n_groups = 15
#  , plot_full_range = T
#)

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.12 [1.08, 1.19]         1.06      0.89
##       barriers_self_cw 1.56 [1.48, 1.66]         1.25      0.64
##       pressure_self_cw 1.06 [1.02, 1.13]         1.03      0.95
##        pushing_self_cw 1.13 [1.09, 1.19]         1.06      0.89
##   facilitators_self_cw 1.81 [1.71, 1.93]         1.35      0.55
##  persuasion_partner_cw 1.28 [1.22, 1.35]         1.13      0.78
##    pressure_partner_cw 1.06 [1.02, 1.13]         1.03      0.95
##     pushing_partner_cw 1.24 [1.19, 1.31]         1.11      0.81
##     persuasion_self_cb 3.74 [3.49, 4.03]         1.93      0.27
##  persuasion_partner_cb 4.03 [3.75, 4.34]         2.01      0.25
##       pressure_self_cb 2.27 [2.13, 2.43]         1.51      0.44
##    pressure_partner_cb 2.22 [2.09, 2.38]         1.49      0.45
##        pushing_self_cb 2.36 [2.21, 2.52]         1.54      0.42
##     pushing_partner_cb 2.43 [2.28, 2.60]         1.56      0.41
##       barriers_self_cb 1.26 [1.20, 1.33]         1.12      0.80
##   facilitators_self_cb 1.14 [1.10, 1.21]         1.07      0.88
##                    day 1.07 [1.04, 1.14]         1.04      0.93
##  Tolerance 95% CI
##      [0.84, 0.92]
##      [0.60, 0.67]
##      [0.89, 0.98]
##      [0.84, 0.92]
##      [0.52, 0.59]
##      [0.74, 0.82]
##      [0.89, 0.98]
##      [0.76, 0.84]
##      [0.25, 0.29]
##      [0.23, 0.27]
##      [0.41, 0.47]
##      [0.42, 0.48]
##      [0.40, 0.45]
##      [0.38, 0.44]
##      [0.75, 0.83]
##      [0.83, 0.91]
##      [0.88, 0.96]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         31%
##          beta         19%
##        cauchy         19%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         75%
##  beta-binomial         16%
##       binomial          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 115 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 486, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.000000000 0.007561553
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summary_is_reactance <- summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.40 0.23 [0.12, 1.25] 0.944 [0.83, 1.20] 0.072 0.353 Weak Evidence for Null 1.000 19541 24420
Within-Person Effects
Daily individual’s experienced persuasion 0.60* 0.14 [0.36, 0.93] 0.987 [0.83, 1.20] 0.066 1.694 Weak Evidence 1.000 18821 21645
Daily partner’s experienced persuasion 1.35 0.36 [0.82, 2.51] 0.877 [0.83, 1.20] 0.295 0.225 Moderate Evidence for Null 1.000 15298 17996
Daily individual’s experienced pressure 6.78* 6.14 [1.40, 63.18] 0.991 [0.83, 1.20] 0.011 3.813 Moderate Evidence 1.000 13931 18160
Daily partner’s experienced pressure 1.51 1.23 [0.29, 13.36] 0.706 [0.83, 1.20] 0.162 0.205 Moderate Evidence for Null 1.000 16765 18592
Daily individual’s experienced pushing 1.59 0.41 [0.99, 2.98] 0.973 [0.83, 1.20] 0.118 0.678 Weak Evidence for Null 1.000 17871 19024
Daily partner’s experienced pushing 0.87 0.28 [0.45, 1.73] 0.670 [0.83, 1.20] 0.391 0.137 Moderate Evidence for Null 1.000 21742 23601
Day 1.65 1.26 [0.37, 7.62] 0.745 [0.83, 1.20] 0.153 0.103 Moderate Evidence for Null 1.000 26409 28870
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.02 0.56 [0.34, 3.35] 0.518 [0.83, 1.20] 0.259 0.096 Strong Evidence for Null 1.000 25169 27753
Daily facilitators 1.08 0.64 [0.30, 3.40] 0.549 [0.83, 1.20] 0.237 0.115 Moderate Evidence for Null 1.000 20365 24604
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 4.44 6.70 [0.24, 116.32] 0.845 [0.83, 1.20] 0.062 0.256 Moderate Evidence for Null 1.000 12271 19664
Mean partner’s experienced persuasion 8.48 13.61 [0.43, 277.62] 0.918 [0.83, 1.20] 0.037 0.371 Weak Evidence for Null 1.000 11662 18399
Mean individual’s experienced pressure 59.85 151.80 [0.54, 15556.44] 0.955 [0.83, 1.20] 0.015 0.708 Weak Evidence for Null 1.000 14113 19622
Mean partner’s experienced pressure 0.73 2.04 [0.00, 144.08] 0.545 [0.83, 1.20] 0.051 0.209 Moderate Evidence for Null 1.000 9900 18663
Mean individual’s experienced pushing 1.59 3.71 [0.02, 178.19] 0.578 [0.83, 1.20] 0.060 0.161 Moderate Evidence for Null 1.000 11104 18225
Mean partner’s experienced pushing 0.01 0.02 [0.00, 1.03] 0.974 [0.83, 1.20] 0.008 1.108 Weak Evidence 1.000 11127 16434
Mean barriers 0.67 0.73 [0.07, 5.64] 0.645 [0.83, 1.20] 0.130 0.129 Moderate Evidence for Null 1.000 18534 26202
Mean facilitators 0.76 0.54 [0.18, 3.15] 0.648 [0.83, 1.20] 0.191 0.174 Moderate Evidence for Null 1.000 14740 22473
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Interactions
persuasion_self_cw:barriers_self_cw 0.76 0.24 [0.38, 1.39] 0.821 [0.83, 1.20] 0.310 0.177 Moderate Evidence for Null 1.000 23607 22535
barriers_self_cw:pressure_self_cw 1.58 1.52 [0.22, 13.05] 0.686 [0.83, 1.20] 0.136 0.175 Moderate Evidence for Null 1.000 23408 28063
barriers_self_cw:pushing_self_cw 0.65 0.25 [0.27, 1.37] 0.880 [0.83, 1.20] 0.194 0.230 Moderate Evidence for Null 1.000 19356 22441
persuasion_self_cw:facilitators_self_cw 0.88 0.25 [0.48, 1.56] 0.679 [0.83, 1.20] 0.436 0.135 Moderate Evidence for Null 1.000 17186 19441
pressure_self_cw:facilitators_self_cw 2.55 3.26 [0.22, 45.77] 0.779 [0.83, 1.20] 0.090 0.321 Weak Evidence for Null 1.000 12990 17446
pushing_self_cw:facilitators_self_cw 0.73 0.31 [0.32, 2.07] 0.758 [0.83, 1.20] 0.232 0.246 Moderate Evidence for Null 1.000 12633 14738
Random Effects
sd(Intercept) 2.64 0.61 [1.64, 4.08] NA NA NA NA NA 1.001 7903 16817
sd(Daily individual’s experienced persuasion) 0.66 0.36 [0.06, 1.45] NA NA NA NA NA 1.000 4439 6807
sd(Daily partner’s experienced persuasion) 0.68 0.35 [0.09, 1.54] NA NA NA NA NA 1.000 7782 7552
sd(Daily individual’s experienced pressure) 2.38 1.02 [0.61, 4.75] NA NA NA NA NA 1.000 8978 7015
sd(Daily partner’s experienced pressure) 1.83 1.17 [0.14, 4.52] NA NA NA NA NA 1.000 11193 12241
sd(Daily individual’s experienced pushing) 0.55 0.40 [0.03, 1.53] NA NA NA NA NA 1.001 6460 10411
sd(Daily partner’s experienced pushing) 0.45 0.37 [0.02, 1.45] NA NA NA NA NA 1.000 11809 15703
sd(barriers_self_cw) 1.13 0.71 [0.08, 2.66] NA NA NA NA NA 1.001 7367 9451
sd(facilitators_self_cw) 1.39 0.60 [0.22, 2.73] NA NA NA NA NA 1.000 7616 7172
sd(persuasion_self_cw:barriers_self_cw) 0.58 0.42 [0.03, 1.61] NA NA NA NA NA 1.001 8107 12898
sd(barriers_self_cw:pressure_self_cw) 0.97 0.85 [0.04, 3.40] NA NA NA NA NA 1.000 14182 15002
sd(barriers_self_cw:pushing_self_cw) 0.84 0.46 [0.08, 1.94] NA NA NA NA NA 1.001 8459 9366
sd(persuasion_self_cw:facilitators_self_cw) 0.33 0.29 [0.02, 1.15] NA NA NA NA NA 1.000 9640 17208
sd(pressure_self_cw:facilitators_self_cw) 2.54 1.21 [0.34, 5.20] NA NA NA NA NA 1.000 10620 8039
sd(pushing_self_cw:facilitators_self_cw) 0.87 0.60 [0.06, 2.36] NA NA NA NA NA 1.001 6284 10446
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.73). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="05_SensitivityBarriersInteraction-ONLY_DAYS_WITH_PLAN-_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw

hypothesis(is_reactance, "exp(pressure_self_cw) > exp(pushing_self_cw)")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (exp(pressure_sel... > 0    11.61        32    -0.03    39.96      18.29
##   Post.Prob Star
## 1      0.95     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = FALSE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
Hurdle Component
Non-Zero Component
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 3.41*** [ 1.79, 6.47] 50.01*** [42.13, 59.09] 123.06*** [108.16, 139.27] 3.79*** [ 3.58, 4.01] 0.40 [0.12, 1.25]
Within-Person Effects
Daily individual’s experienced persuasion 1.87*** [ 1.41, 2.61] 1.00 [ 0.95, 1.06] 1.02 [ 0.98, 1.06] 0.01 [-0.04, 0.06] 0.60* [0.36, 0.93]
Daily partner’s experienced persuasion 1.40** [ 1.12, 1.83] 1.02 [ 0.98, 1.07] 1.03 [ 0.99, 1.07] 0.01 [-0.04, 0.06] 1.35 [0.82, 2.51]
Daily individual’s experienced pressure 0.86 [ 0.39, 2.09] 0.91 [ 0.78, 1.05] 0.97 [ 0.87, 1.07] -0.02 [-0.16, 0.12] 6.78* [1.40, 63.18]
Daily partner’s experienced pressure 1.69 [ 0.85, 5.72] 0.97 [ 0.88, 1.06] 0.97 [ 0.90, 1.05] -0.05 [-0.18, 0.07] 1.51 [0.29, 13.36]
Daily individual’s experienced pushing 1.24 [ 0.77, 2.01] 0.99 [ 0.93, 1.06] 1.02 [ 0.96, 1.08] 0.02 [-0.06, 0.10] 1.59 [0.99, 2.98]
Daily partner’s experienced pushing 1.60** [ 1.13, 2.54] 0.96 [ 0.90, 1.02] 1.00 [ 0.96, 1.05] 0.09* [ 0.02, 0.16] 0.87 [0.45, 1.73]
Day 0.91 [ 0.55, 1.51] 0.95 [ 0.83, 1.09] 0.96 [ 0.87, 1.06] 0.26*** [ 0.11, 0.41] 1.65 [0.37, 7.62]
Own actionplan NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.14 [ 0.64, 2.04] 1.26*** [ 1.13, 1.41] 1.01 [ 0.94, 1.07] -0.12 [-0.25, 0.01] 1.02 [0.34, 3.35]
Daily facilitators 1.17 [ 0.68, 1.96] 1.19** [ 1.06, 1.33] 1.04 [ 0.99, 1.10] 0.11** [ 0.05, 0.17] 1.08 [0.30, 3.40]
Daily wear time NA NA NA NA 1.00*** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.07 [ 0.27, 4.02] 1.06 [ 0.75, 1.48] 0.99 [ 0.69, 1.40] 0.13 [-0.47, 0.72] 4.44 [0.24, 116.32]
Mean partner’s experienced persuasion 1.43 [ 0.38, 5.24] 0.96 [ 0.69, 1.32] 0.86 [ 0.61, 1.23] 0.42 [-0.18, 1.00] 8.48 [0.43, 277.62]
Mean individual’s experienced pressure 0.18 [ 0.03, 1.11] 0.68 [ 0.29, 1.58] 1.12 [ 0.77, 1.64] -0.07 [-0.67, 0.53] 59.85 [0.54, 15556.44]
Mean partner’s experienced pressure 0.19 [ 0.03, 1.17] 0.49 [ 0.21, 1.15] 1.04 [ 0.72, 1.49] -0.32 [-0.91, 0.28] 0.73 [0.00, 144.08]
Mean individual’s experienced pushing 1.32 [ 0.20, 8.70] 1.53 [ 0.84, 2.80] 1.11 [ 0.69, 1.82] 0.23 [-0.61, 1.08] 1.59 [0.02, 178.19]
Mean partner’s experienced pushing 1.45 [ 0.22, 9.51] 1.65 [ 0.91, 3.01] 1.27 [ 0.78, 2.05] -0.06 [-0.90, 0.79] 0.01 [0.00, 1.03]
Mean barriers 0.50** [ 0.30, 0.83] 0.97 [ 0.83, 1.12] 0.80*** [ 0.72, 0.89] -0.50*** [-0.65, -0.34] 0.67 [0.07, 5.64]
Mean facilitators 1.24 [ 0.86, 1.80] 1.08 [ 0.98, 1.19] 0.96 [ 0.89, 1.03] 0.29*** [ 0.18, 0.39] 0.76 [0.18, 3.15]
Mean wear time NA NA NA NA 1.00* [ 1.00, 1.00] NA NA NA NA
Interactions
persuasion_self_cw:barriers_self_cw 1.02 [ 0.65, 1.63] 1.03 [ 0.96, 1.11] 0.97 [ 0.93, 1.02] 0.00 [-0.08, 0.07] 0.76 [0.38, 1.39]
barriers_self_cw:pressure_self_cw 0.34 [ 0.09, 1.21] 0.95 [ 0.75, 1.20] 0.96 [ 0.81, 1.14] -0.14 [-0.47, 0.20] 1.58 [0.22, 13.05]
barriers_self_cw:pushing_self_cw 2.50* [ 1.14, 6.03] 0.94 [ 0.85, 1.05] 1.02 [ 0.95, 1.09] 0.08 [-0.02, 0.18] 0.65 [0.27, 1.37]
persuasion_self_cw:facilitators_self_cw 1.29 [ 0.87, 2.10] 1.02 [ 0.96, 1.09] 0.99 [ 0.95, 1.04] 0.00 [-0.06, 0.06] 0.88 [0.48, 1.56]
pressure_self_cw:facilitators_self_cw 0.88 [ 0.18, 4.87] 0.94 [ 0.73, 1.18] 0.92 [ 0.69, 1.14] 0.00 [-0.24, 0.31] 2.55 [0.22, 45.77]
pushing_self_cw:facilitators_self_cw 1.54 [ 0.77, 3.92] 0.95 [ 0.88, 1.04] 1.03 [ 0.96, 1.09] 0.02 [-0.08, 0.13] 0.73 [0.32, 2.07]
Random Effects
sd(Intercept) 1.64 [1.23, 2.23] 0.32 [0.25, 0.43] 0.34 [0.26, 0.45] 0.62 [0.49, 0.81] 2.64 [1.64, 4.08]
sd(Daily individual’s experienced persuasion) 0.21 [0.01, 0.61] 0.08 [0.03, 0.14] 0.05 [0.01, 0.09] 0.04 [0.00, 0.11] 0.66 [0.06, 1.45]
sd(Daily partner’s experienced persuasion) 0.30 [0.03, 0.66] 0.06 [0.01, 0.11] 0.05 [0.01, 0.10] 0.04 [0.00, 0.11] 0.68 [0.09, 1.54]
sd(Daily individual’s experienced pressure) 0.54 [0.02, 2.15] 0.08 [0.00, 0.33] 0.05 [0.00, 0.20] 0.07 [0.00, 0.25] 2.38 [0.61, 4.75]
sd(Daily partner’s experienced pressure) 0.77 [0.04, 2.37] 0.05 [0.00, 0.18] 0.03 [0.00, 0.12] 0.09 [0.00, 0.29] 1.83 [0.14, 4.52]
sd(Daily individual’s experienced pushing) 0.70 [0.07, 1.53] 0.05 [0.00, 0.14] 0.08 [0.01, 0.15] 0.09 [0.01, 0.20] 0.55 [0.03, 1.53]
sd(Daily partner’s experienced pushing) 0.44 [0.03, 1.23] 0.06 [0.00, 0.14] 0.04 [0.00, 0.11] 0.09 [0.01, 0.19] 0.45 [0.02, 1.45]
sd(barriers_self_cw) 1.51 [1.10, 2.08] 0.23 [0.14, 0.35] 0.13 [0.06, 0.20] 0.32 [0.23, 0.43] 1.13 [0.08, 2.66]
sd(facilitators_self_cw) 1.30 [0.84, 1.95] 0.27 [0.19, 0.39] 0.07 [0.01, 0.15] 0.05 [0.00, 0.16] 1.39 [0.22, 2.73]
sd(persuasion_self_cw:barriers_self_cw) 0.58 [0.05, 1.26] 0.07 [0.00, 0.18] 0.03 [0.00, 0.09] 0.07 [0.00, 0.18] 0.58 [0.03, 1.61]
sd(barriers_self_cw:pressure_self_cw) 0.67 [0.03, 2.58] 0.11 [0.01, 0.42] 0.07 [0.00, 0.31] 0.41 [0.08, 0.85] 0.97 [0.04, 3.40]
sd(barriers_self_cw:pushing_self_cw) 1.26 [0.38, 2.39] 0.15 [0.05, 0.26] 0.03 [0.00, 0.10] 0.05 [0.00, 0.18] 0.84 [0.08, 1.94]
sd(persuasion_self_cw:facilitators_self_cw) 0.59 [0.11, 1.23] 0.07 [0.00, 0.16] 0.03 [0.00, 0.09] 0.04 [0.00, 0.13] 0.33 [0.02, 1.15]
sd(pressure_self_cw:facilitators_self_cw) 1.60 [0.10, 4.28] 0.12 [0.01, 0.51] 0.24 [0.02, 0.62] 0.16 [0.01, 0.55] 2.54 [0.34, 5.20]
sd(pushing_self_cw:facilitators_self_cw) 0.94 [0.08, 2.45] 0.05 [0.00, 0.17] 0.04 [0.00, 0.14] 0.11 [0.01, 0.26] 0.87 [0.06, 2.36]
Additional Parameters
sigma NA NA 0.62 [0.59, 0.65] 0.53 [0.51, 0.55] 0.84 [0.82, 0.88] NA NA

Plot Interaction effect.

  1. Probability of being active (crossing the hurdle)
library(ggplot2)
library(plotly)


p <- conditional_effects(
  pa_sub, 
  dpar = 'hu',
  effects = 'pushing_self_cw:barriers_self_cw',
  re_formula = NA, 
  method = 'posterior_linpred',
  plot = FALSE
)

p2 <- p

p2[[1]]$estimate__ <- 1 - p[[1]]$estimate__
p2[[1]]$lower__ <- 1 - p[[1]]$upper__
p2[[1]]$upper__ <- 1 - p[[1]]$lower__

plots <- plot(p2, theme = theme_minimal(), plot = FALSE) 

plots[[1]] + 
  scale_color_manual(values = c("steelblue", "forestgreen", "coral")) +
  scale_fill_manual(values = c("steelblue", "forestgreen", "coral"))

Test slope of pushing at high levels of barriers

library(tidybayes)
library(dplyr)
library(ggplot2)

plot_slope_by_moderator <- function(
  mod,            
  df,             
  moderator,      
  predictor_par,  
  interaction_par,
  exponentiate = FALSE,   
  n_points = 1000,
  fill_colors = c("Significant" = "steelblue",
                  "Non-significant" = "coral"),
  line_colors = fill_colors,  # Reuse fill colors by default
  vjust_val = -1.5
) {
  
  # Generate a sequence of moderator values
  moderator_seq <- seq(
    from = min(df[[moderator]], na.rm = TRUE),
    to = max(df[[moderator]], na.rm = TRUE),
    length.out = n_points
  )
  
  # Extract posterior draws for the parameters of interest
  library(rlang)
  predictor_sym <- rlang::sym(predictor_par)
  interaction_sym <- rlang::sym(interaction_par)

  post_draws <- mod %>%
    spread_draws(!!predictor_sym, !!interaction_sym)
  
  # For each moderator value, compute the slope for each posterior draw
  slope_list <- lapply(moderator_seq, function(mod_val) {
    tibble(
      moderator_value = mod_val,
      slope = - (post_draws[[predictor_par]] + mod_val * post_draws[[interaction_par]])
    )
  })
  
  # Combine the list into a single data frame
  slope_df <- bind_rows(slope_list)
  
  # Summarize the posterior for each moderator value (using a 95% credible interval)
  slope_summary <- slope_df %>%
    group_by(moderator_value) %>%
    summarize(
      slope_median = median(slope),
      slope_lower  = quantile(slope, probs = 0.025),
      slope_upper  = quantile(slope, probs = 0.975)
    ) %>%
    # Create significance flag: "sig" if CI excludes 0, "ns" otherwise
    mutate(sig_flag = if_else(slope_lower > 0 | slope_upper < 0, "sig", "ns")) %>%
    arrange(moderator_value)
  
  # If exponentiate == TRUE, transform slope columns
  if (exponentiate) {
    slope_summary <- slope_summary %>%
      mutate(
        slope_median_exp = exp(slope_median),
        slope_lower_exp  = exp(slope_lower),
        slope_upper_exp  = exp(slope_upper)
      )
  }
  
  # Identify boundary rows and duplicate them for each category
  slope_summary <- slope_summary %>%
    mutate(
      next_sig_flag = lead(sig_flag),
      is_boundary   = sig_flag != next_sig_flag & !is.na(next_sig_flag)
    )
  
  boundary_rows <- slope_summary %>%
    filter(is_boundary)
  
  boundary_rows_next <- boundary_rows
  boundary_rows_next$sig_flag <- boundary_rows_next$next_sig_flag
  
  slope_summary_extended <- bind_rows(slope_summary, boundary_rows_next)
  
  # Split into "sig" and "ns"
  slope_summary_sig <- slope_summary_extended %>% filter(sig_flag == "sig")
  slope_summary_ns  <- slope_summary_extended %>% filter(sig_flag == "ns")
  
  # Build the base plot:
  #    We'll use either slope_{median, lower, upper} or slope_{median_exp, lower_exp, upper_exp} 
  #    depending on exponentiate.
  if (exponentiate) {
    median_col <- "slope_median_exp"
    lower_col <- "slope_lower_exp"
    upper_col <- "slope_upper_exp"
    y_lab <- "Odds Ratio (Exponentiated Slope)"
    ref_line <- 1  # reference line for ORs is 1
  } else {
    median_col <- "slope_median"
    lower_col <- "slope_lower"
    upper_col <- "slope_upper"
    y_lab <- "Slope (Log-Odds Scale)"
    ref_line <- 0
  }
  
  p <- ggplot() +
    # Non-significant region
    geom_ribbon(
      data = slope_summary_ns,
      aes_string(
        x = "moderator_value",
        ymin = lower_col,
        ymax = upper_col,
        fill = '"Non-significant"'
      ),
      alpha = 0.2
    ) +
    geom_line(
      data = slope_summary_ns,
      aes_string(
        x = "moderator_value",
        y = median_col,
        color = '"Non-significant"'
      ),
      size = 1
    ) +
    # Significant region
    geom_ribbon(
      data = slope_summary_sig,
      aes_string(
        x = "moderator_value",
        ymin = lower_col,
        ymax = upper_col,
        fill = '"Significant"'
      ),
      alpha = 0.2
    ) +
    geom_line(
      data = slope_summary_sig,
      aes_string(
        x = "moderator_value",
        y = median_col,
        color = '"Significant"'
      ),
      size = 1
    ) +
    # Reference line (0 for log-odds, 1 for odds-ratios)
    geom_hline(
      yintercept = ref_line, 
      linetype = "dashed"
    ) +
    # Manual scales for fill & color
    scale_fill_manual(
      values = fill_colors
    ) +
    scale_color_manual(
      values = line_colors
    ) +
    labs(
      x = moderator,
      y = y_lab,
      fill = "Significance",
      color = "Significance"
    ) +
    theme_minimal()
  
  # Add boundary lines (where significance changes) & labels
  p2 <- p +
    geom_vline(
      data = boundary_rows,
      aes(xintercept = moderator_value),
      color = fill_colors["Significant"],
      alpha = 0.35,
      size = 0.7,
      linetype = "dashed"
    ) +
    geom_text(
      data = boundary_rows,
      aes(x = moderator_value, y = -Inf, label = round(moderator_value, 3)),
      vjust = vjust_val,
      size  = 3.7
    )
  
  return(p2)
}


p_log_odds <- plot_slope_by_moderator(
  mod = pa_sub, 
  df = df_double, 
  moderator = "barriers_self_cw", 
  predictor_par = "b_hu_pushing_self_cw", 
  interaction_par = "b_hu_barriers_self_cw:pushing_self_cw",
  exponentiate = FALSE
)

p_log_odds

p_odds_ratio <- plot_slope_by_moderator(
  mod = pa_sub, 
  df = df_double, 
  moderator = "barriers_self_cw", 
  predictor_par = "b_hu_pushing_self_cw", 
  interaction_par = "b_hu_barriers_self_cw:pushing_self_cw",
  exponentiate = TRUE,
  vjust_val = -12,
) + coord_cartesian(ylim = c(1, 18))

p_odds_ratio

# Export
ggsave(
  filename = file.path("Output", "InteractPlot_log_odds.png"), 
  plot = p_log_odds,
  width = 8,
  height = 5,
  dpi = 300,
  bg = 'white'
)

ggsave(
  filename = file.path("Output", "InteractPlot_odds_ratio.png"), 
  plot = p_odds_ratio, 
  width = 8, 
  height = 5, 
  dpi = 300,
  bg = 'white'
)
  1. Expected MVPA when the hurdle was crossed (interaction NOT significant)
p3 <- conditional_effects(
  pa_sub, 
  dpar = 'mu',
  effects = 'pushing_self_cw:barriers_self_cw',
  #re_formula = NULL, 
  plot = FALSE
)

p3[[1]]$estimate__ <- exp(p3[[1]]$estimate__)
p3[[1]]$lower__ <- exp(p3[[1]]$lower__)
p3[[1]]$upper__ <- exp(p3[[1]]$upper__)

plot(p3, ask = FALSE)

  1. Expected physical activity in Minutes (when both components combined)
# Expected physical activity
p4 <- conditional_effects(
  pa_sub, 
  effects = 'pushing_self_cw:barriers_self_cw',
  #re_formula = NULL,
  plot = FALSE
)

plot(p4, ask = FALSE)

report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()
  • rmarkdown (version 2.29; Allaire J et al., 2024)
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • rlang (version 1.1.4; Henry L, Wickham H, 2024)
  • tidybayes (version 3.0.7; Kay M, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.2; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • plotly (version 4.10.4; Sievert C, 2020)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()